# Code to create Figure A1 in "Pandemics and Cities: Evidence from the Black Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 5-14-20
# Last updated: 9-5-23


library(sf)
library(raster)
library(tidyverse)
library(tmap)
library(haven)

# Need the following features in the map:

# (1) Modern Country Borders
# (2) Cities: ASTRAHAN, FEODOSIJA, MESSINA, GENOVA, IZMAIL, WIEN, MOSKVA, NOVGOROD, VISBY, Prague, Leipzig, Halych
# (3) Rivers: Danube, Dniester, the river from Astrakhan to Moscow, moscow to Novgord
# (4) Medieval Trade Routes: the medieval land route (Via Regia), 
# (5) seas

# coordinates of the cities not in Bosker (from wikipedia)
# c("ASTRAHAN",51.166667, 71.433333)
# c("FEODOSIJA",45.048889, 35.379167)
# c("VINICA",45.351667, 28.836389)
# c("MOSKVA",55.755833, 37.617222)
# c("NOVGOROD",58.55, 31.266667)
# c("VISBY",57.634722, 18.299167)

# set working directory
setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A1")

# bring in the cities
cities_raw <- read_dta("All_Cities.dta")
# keep the variables you want
cities_raw <- cities_raw %>% select(city_jjk, longitude, latitude)
# keep the cities you want
target <- c("WIEN", "PRAHA", "LEIPZIG")
cities_raw <- cities_raw %>% filter(city_jjk %in% target)
# add the cities not contained in Bosker that you acquired coordinates for form Bosker
cities_raw <- cities_raw %>%
  rbind(c("ASTRAHAN (1345)",48.05,46.35)) %>%
  rbind(c("KAFFA (1346)",35.379167,45.048889)) %>%
  rbind(c("VINICA",28.836389,45.351667)) %>%
  rbind(c("MOSKVA",37.617222,55.755833)) %>%
  rbind(c("NOVGOROD",31.266667,58.55)) %>%
  rbind(c("VISBY",18.299167,57.634722)) %>%
  rbind(c("HALYCH",24.728611,49.124722)) %>%
  rbind(c("GENOVA",8.946256,44.4056499)) %>%
  rbind(c("MESSINA (1347)",15.5540152,38.1938137))
# clean up the cities data
glimpse(cities_raw)
cities_raw$longitude <- as.numeric(as.character(cities_raw$longitude))
cities_raw$latitude <- as.numeric(as.character(cities_raw$latitude))
glimpse(cities_raw)

# spatialize the cities data
cities <- st_as_sf(cities_raw, coords=c("longitude", "latitude"), crs=4326)

# take a look at the data
cities
glimpse(cities)
head(cities)
plot(cities[2], col="red", pch=16, reset=F)

# Make a bounding box of the cities so you can create one a little bigger
bbox_cities <- st_bbox(cities, crs=4326)
bbox_cities
xrange <- bbox_cities$xmax - bbox_cities$xmin # range of x values
yrange <- bbox_cities$ymax - bbox_cities$ymin # range of y values
bbox_cities[1] <- bbox_cities[1] - (0.20 * xrange) # xmin - left
bbox_cities[3] <- bbox_cities[3] + (0.20 * xrange) # xmax - right
bbox_cities[2] <- bbox_cities[2] - (0.20 * yrange) # ymin - bottom
bbox_cities[4] <- bbox_cities[4] + (0.20 * yrange) # ymax - top
# crop the cities using your bounding box
cities <- cities %>%
  st_crop(bbox_cities)

# Bring in medieval roads shape file
med_roads <- st_read("Europe_Medieval_Trade_Routes/MedRdsProj.shp") %>%
  st_transform(4326)
# take a look at the data
glimpse(med_roads)
head(med_roads)
# get rid of everything except the gis field
med_roads <- med_roads %>% dplyr::select()
# crop the roads
med_roads <- med_roads %>%
  st_crop(bbox_cities)
glimpse(med_roads)
# add Via Reggia
crs = st_crs(4326)
s1 <- rbind(c(24.728611,49.124722),c(14.4378005, 50.0755381))
ls <- st_linestring(s1)
b = st_sf(geometry = st_sfc(ls), crs = crs)
med_roads <- med_roads %>%
  rbind(b)
# plot the roads
plot(med_roads[1], col="brown", reset=F, add=T)

# Bring in rivers shape file: round 1
rivers <- st_read("Rivers/MajorRivers.shp") %>%
  st_transform(4326)
glimpse(rivers)
rivers <- rivers %>% dplyr::select(NAME)
rivers <- rivers %>%
  st_crop(bbox_cities)
glimpse(rivers)
plot(rivers[2], col="blue", reset=F, add=T)

# Bring in rivers shape file: round 2
rivers2 <- st_read("Rivers/Rivers_conic.shp") %>%
  st_transform(4326)
glimpse(rivers2)
rivers2 <- rivers2 %>% dplyr::select(name)
rivers2 <- rivers2 %>%
  st_crop(bbox_cities)
glimpse(rivers2)
target <- c("Dniester", "Volga", "Oka")
rivers2 <- rivers2 %>% filter(name %in% target)
plot(rivers2[2], col="blue", reset=F, add=T)

# Bring in the modern country borders
modern_countries <- st_read("ModernCountries/Longitude_Graticules_and_World_Countries_Boundaries-shp/99bfd9e7-bb42-4728-87b5-07f8c8ac631c2020328-1-1vef4ev.lu5nk.shp") %>%
  st_transform(4326)
# Clip the countries using bounding box based on the cities data
modern_countries <- modern_countries %>%
  st_crop(bbox_cities)
glimpse(modern_countries)
head(modern_countries)
modern_countries <- modern_countries %>% dplyr::select(CNTRY_NAME)
glimpse(modern_countries)
plot(modern_countries[2], col="NA", reset=F, add=T)

# Bring in seas shape file
seas <- st_read("Seas/ne_10m_ocean.shp") %>%
  st_transform(4326)
glimpse(seas)
# there is a ring intersection error in the seas .shp file. So need to fix the geometry.
sf::sf_use_s2(FALSE) # this does work
seas_crop <- st_crop(seas, bbox_cities)
sf::sf_use_s2(TRUE) # putting settings back to their default



# make figure A1
state_map <- tm_shape(modern_countries) + tm_borders(col = "grey") 

rivers_map <- tm_shape(rivers) +
  tm_lines(col="blue", lwd=3) +
  tm_add_legend("line", col = "blue", title = "Major Rivers", lwd=3)

rivers_map2 <- tm_shape(rivers2) +
  tm_lines(col="blue", lwd=3)

med_map <- tm_shape(med_roads) +
  tm_lines(col="black") +
  tm_add_legend("line", col = "black", title = "Trade Routes")
 
seas_map <- tm_shape(seas_crop) +
  tm_polygons(col="grey")

cities_map <- tm_shape(cities) +
  tm_text("city_jjk", size=0.75,ymod=1,xmod=.75, fontface="bold") +
  tm_dots(col="red",size=.40) +
  tm_add_legend("symbol", col = "red", title = "Trade Cities")

full_map <- rivers_map + rivers_map2 + med_map + cities_map + seas_map
full_map <- full_map + tm_layout(legend.position = c(.77,.70),legend.frame=T)
full_map

# save the map (note: the orange arrows were inputted by hand)
tmap_save(full_map, filename="Trade_Figure1.png",
          height=8.5, width=11, units="in", dpi=300)

#











# End Code